The Role of Substrate Corrugation in Helium Monolayer Solidification 
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We investigate the first layer of helium adsorbed on graphite with path-integral Monte Carlo, 
examining the role of substrate corrugations on the phase diagram. When no corrugations are 
present, the equilibrium state of the system is a liquid phase, with solidification occurring only 
under compression but before layer promotion. We determine the solid- liquid coexistence region and 
compare our results to recent Green's function Monte Carlo calculations on the same system. When 
substrate corrugations are included, we find that the equilibrium phase is the y/3x y/s commensurate 
solid phase that is well known from experiment. The melting behavior, heat capacity, and single 
particle binding energy are determined and compared to experiment. We further find that for 
densities below the commensurate coverage, the low temperature phase of the system consists of 
solid clusters in coexistence with coalesced vacancies. We find no first layer liquid phase and so no 
superfluidity in this layer, in contrast to some rather recent suggestions. 



PACS numbers 67.70. +n, 67.40 Kh 
I. INTRODUCTION 

Quantum films such as helium adsorbed on graphite 
surface are characterized by a rich phase diagram as one 
varies the number of the deposited helium atoms. As a 
function of the helium coverage, the helium film grows 
in atomically thin layerstrEl and at least seven such lay- 
ers may be clearly distinguished and studied on a well- 
prepared substrate. 

During the iaat several years, extensive heat capacity 
measurementsaH of the first six layers have been per- 
formed, and superfluidity in the higher layers has been 
detected by both torsional oscillatorLrQ and third sound 
measurements. More specifically, the first layer of helium 
adsorbed on graphite has been the laboratory for study 
of a variety of phenomena. In this layer a V3 x y/3 com- 
mensurate solid phase forms in which one-third , 
available substrate adsorption sites are occupicd.l 
For densities above this commensurate density, the mono- 
layer is characterized by a region of domain wall phases 
and at even higher monolayer densitieSpit forms an in- 
commensurate triangular solid phase.EiTEj The phase di- 
agram below the commensurate density at low tempera- 
tures is not well understood and there are two competing 
scenarios. In one scenario this— region corresponds to a 
solid with clustered vacancies £.3 According to this pic- 
ture, since the commensurate solid phase is in th£ same 
universality class as the three-state Potts model,Ej~t3 at 
lower densities the film should consist of a commensurate 
solid with vacancies. When the temperature is raised, 
the solid melts continuously. Lowering the temperature 
causes the vacancies to coalesce, which is a first order 
transition. The difference between the temperatures of 
these two transitions becomes smaller as the density is 
lowered until they meet at a tricritical point. 

The second scenario for the low density region of the 
phase diagram of the sub- monolayer was suggested more 
recently by Grey wall and Buschcl (GB). GB point out 



that their measured heat capacity is not linear in density 
for the entire region below the commensurate density, 
as it must be for phase coexistence. Thus, they pro- 
pose that a self-bound liquid phase occurs at about 0.04 
atom/ A 2 . This conclusion is supported by 2ELvai : iAtional 
calculations .E3 However, direct measurementso'BO detect 
no superfluidity, possibly because of poor substrate con- 
nectivity. | — | 

In a previous publications, we have used the path- 
integral Monte Carlo method and realistic helium-helium 
and helium graphite interactions to study the monolayer 
at or below the commensurate density. We found that 
the presence of corrugations, which causes the commen- 
surate solid at 1/3 coverage, creates solid commensurate 
clusters at densities below the commensurate density. We 
found that it is the melting of these monolayer clusters 
which gives rise to a specific heat maximum which was 
incorrectly interpretecH as onset of monolayer superflu- 
idity. In this paper we examine in detail the role of sub- 
strate corrugation on the first layer phase diagram using 
the path-integral Monte Carlo method. First, we present 
results for helium adsorbed on a smooth substrate. The 
helium-graphite interaction is based-pn the laterally av- 
eraged potential of Carlos and Cole£3 This layer exhibits 
liquid and solid phases, and we calculate the phase dia- 
gram at low temperatures. These results are compared, 
with recent Green's function Monte Carlo calculations .til 
The second layer promotion density is also determined. 
Next, calculations including substrate corrugations are 
discussed, and direct comparisons of the melting behav- 
ior, specific heat, and single particle binding energy that 
we obtain are made with experiment. Finally, having 
verified our method by the above comparisons, we exam- 
ine the low density, low temperature phase of the helium 
monolayer. We determine that for all coverages below 
the commensurate density, the system consists of solid 
clusters in phase coexistence with coalesced vacancies. 
No liquid phase occurs and so there is no possibility for 
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first layer superfluidity. 



II. SIMULATION METHOD 

Our study employs a path-integral Monte Carlo 
(PIMC) method for simulating bose systems. Details on 
the application of this method j— to bulk helium can be 
found in the review by Ceperley,E3 and our modifications 
for the simulation of adsorbed helium films may be found 
in a previous publication.E£l We briefly review the proce- 
dure now in order to explain the calculations presented 
in this paper. 

PIMC evaluates properties of an A^-body quantum sys- 
tem at the inverse temperature r by sampling the par- 
tition function Z . Z is expanded as a path-integral by 
inserting M intermediate configurations: 
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Xp(R lf R 2 ; r)p(R 2 , R 3 ; r) . . . p(R M , PR l5 r), (1) 

where p is the density matrix, R^ is a configuration of N 
particles, and r = (3/M . By taking M large enough, the 
density matrices at r can be accurately approximated. 
The particular PIMC method that we employ ergodi- 
cally samples both particle positions and particle permu- 
tations. 

The key element needed in the first layer simulation is 
an accurate approximation for the high temperature den- 
sity matrices. For bulk helium and helium on a smooth 
substrate, the starting approximation for the density ma- 
trix can be taken with a temperature as low as 40 K. In 
this paper, results for the first layer on a smooth sub- 
strate were obtained using the same starting approxima- 
tion for the density matrix that— .was used in previous 
calculations for the second layer .l3 |— , 

However, as discussed in our previous publication) 2 ^ a 
helium-graphite interaction that includes substrate cor- 
rugations makes the starting approximation used for 
smooth substrates impractical, and so we must use a 
simpler form. For these calculations, we instead use the 
high temperature expansion of the density matrix. At 
sufficiently small values of the inverse temperature, the 
density matrix is given by the Trotter approximation: 



exp(— tH) exp(— tT) exp(— tV), 



(2) 



where H, T and V are the Hamiltonian and the kinetic 
and potential energy operators, respectively, and r is the 
inverse temperature. This allows us to approximate the 
density matrices at r as 

p(R,R';r) cx exp[~7r(R - R +1 )/A T 

-T(V(R i ) + V(R i+1 ))/2], (3) 

where R and R' are the position vectors for N particle 
configurations, and At is the thermal wavelength. This 



is sometimes referred to as the semiclassical approxima- 
tion, and is accurate for sufficiently high starting tem- 
peratures. Averaging over the potential energy terms is 
referred to as the endpoint approximation. The potential 
energy term V(Ri) is the sum of all helium- helium and 
helium-graphite interactions. The exponent ,of Eq. |^ is 
the first term in an expansion in powers of tE3 

In the semiclassical calculations, we have used 200 K 
as the starting temperature, meaning that 200 inverse- 
temperature slices are required to reach 1 K. We have 
verified that this temperature is sufficiently high for the 
approximation to be accurate by comparing the energy 
calculated at 4 K using 200 K and 320 K as the starting 
temperatures. The energy values obtained agreed within 
error bars. We have also adopted a three-level bisection, 
I = 3, for sampling the positions. We verified that this 
level gave a lower energy at 4 K than calculations using 
I = 2 and I = 4. The acceptance rate using I = 3 is 
approximately 40%. We have further verified that the 
density matrix is well approximated at 200 K by Eq. || 
by including the next term in the series expansion in 
powers of r. Energy values calculated with this starting 
approximation agreed within error bars with those that 
used Eq. ||. 

We did find that higher order terms in r must be in- 
cluded in the energy estimator. The values we report in 
Sec. IV are obtained using 
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<E>= < T > + < V > +T 2 h 2 /(8m) < (W) 2 > (4) 

2t 

for the energy expectation value, where N is the number 
of particles, < T > and < V > are the kinetic and po- 
tential energy expectation values, and m is the mass of 
the helium atom. Note that in PIMC the effective action, 
defined as the natural log of the density matrix, not the 
energy estimator, is used to choose between configura- 
tions in the Monte Carlo procedure. 

Once the machinery for the Monte Carlo method is in 
place, helium-helium and helium-substrate potentials are 
required for input. We use the AzizE3 potential for the 
helium-helium interaction. For the helium-graphite inter- 
action, we have adopted the anisotropic, Lennard- Jones 
potential proposed by Carlos and Cole.E 3 ] This potential 
can be expressed as a Fourier series in the reciprocal lat- 
tice vectors, G, of the graphite substrate. In cylindrical 
coordinates (/3,z) the expansion is 



V(r) = V (z) + }2 Vg(z) exp(zG ■ p), 

G 



(5) 



where Vq (z) is the laterally averaged potential, and Vg (z) 
gives the corrugation strength. Tha-^nathematical forms 
for these terms are given elsewhere J23 and the series con- 
verges rapidly. For the smooth substrate, we use only 
Vo(z). For calculations that included corrugations, we 
kept the six lowest, equivalent values of G in the expan- 
sion. 
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The limitation of our method is thus related to the ac- 
curacy of the potentials available to us. It is possible, 
for instance, that the substrate may substantially medi- 
ate the helium-helium iateraction.cj This is the so-called 
McLachlan interaction^. We have performed calcula- 
tions both with and without this term. Another possible 
concern is the helium-graphite potential, Eq. ^, which 
may overestimate the corrugation strength. Lowering 
the corrugation height will effect the properties of the 
first layer, with the favored phase becoming liquid in- 
stead of solid at sufficiently small corrugations. We have 
examined this effect by repeating some of the calculations 
with the corrugation strength set to 50 % and to 75 % of 
the value obtained from the Carlos-Cole model. Details 
of these tests of the limits of our interaction model are 
given below. 

Particle permutations were also included in the cal- 
culations and were observed in the film on the smooth 
substrate at low densities. However, we have found that 
permutations do not play a role in the first layer on the 
corrugated substrate. We have allowed for permutations 
at intermediate densities for temperatures as low as 0.571 
K with the level of the path bisectioning taken as large as 
5 (32 beads updated in one move) but never observed any 
particle exchanges. This has been checked in calculations 
both with and without the McLachlan interaction. 



III. RESULTS FOR SMOOTH GRAPHITE 
SURFACE 

A. Energy Calculations 

In this section we report simulation results for the first 
layer obtained by using only the laterally averaged por- 
tion Vq(z) of the helium-graphite potential. We ignore 
substrate corrugations. This system has recently been 
studied jsdth the Green's Function Monte Carlo (GFMC) 
method,Eil so we have the opportunity to verify that the 
full implementation of our method is in agreement with 
the results of these complementary calculations. The 
GFMC calculations employ the same helium-graphite po- 
tential as we do.— but use an older form of the helium- 
helium potentially The older form was used since the au- 
thors wanted to make direct comparisons with previous, 
work on two-dimensional helium using this potential.^ 
The newer potential that we use for the helium-helium 
interaction is somewhat more attractive (the well depth 
is about 0.1 K greater), so we expect the energy per par- 
ticle to be somewhat lower near the equilibrium density. 
This has been observed in recent zero-temperature cal- 
culations for two dimensional liquid helium that employ 
the newer potentials 

Guided by previous simulations, we expect the helium 
film to have a self-bound liquid phase and to solidify at 
high densities, before promotion to the next layer occurs. 
Calculations at low and intermediate densities are per- 



formed using a square simulation cell, while solid phase 
calculations use a rectangular simulation cell that accom- 
modates the periodic structure of the triangular solid. It 
is not necessary in PIMC to employ different forms for 
the density matrix for the liquid and solid phases. Liq- 
uid calculations employed 16 particles, while solid phase 
calculations used 30 particles. Both sets of calculations 
were performed at 400 mK. Calculations at 500 mK are 
in agreement with these results, indicating that we have 
converged to the zero-temperature limit. We verified that 
there were no finite-size errors in the liquid phase by re- 
peating some of the calculations with 32 particles. The 
energy values at the two temperatures were in agreement. 
Finite size errors were found to be negligible for the solid 
phase also, since the energy value calculated at 0.0689 
atom/ A 2 using the rectangular simulation cell (30 par- 
ticles) agreed with the value calculated using the square 
simulation cell (16 particles). 

These calculations are somewhat different from those 
we discuss in the other section of this paper. The poten- 
tial between the active helium atoms and the underlying 
substrate is featureless in the plane of the substrate, so 
the size of the simulation cell may be varied continu- 
ously. This allows us to keep the number of particles 
constant for each phase. In contrast, other calculations 
we have performed (including those with substrate cor- 
rugations taken into account) were with a varying num- 
ber of particles and a constant simulation cell size. The 
two methods lead to somewhat different forms for the 
Maxwell construction. In the present case, liquid-solid 
coexistence will be characterized by a linear region in 
the dependence of the energy per particle on the atomic 
area (inverse volume). 

Figure [l] displays our results for the energy of the first 
layer liquid. Also shown are the results obtained from 
GFMC. For all points except the near equilibrium, the 
two calculations agree within error bars. Near the energy 
minimum, there is some disagreement, but even here the 
computed values differ only by about 0.6%. This is per- 
haps attributable to the different helium-helium poten- 
tials used in the two calculations. Figure || shows our 
results and the GFMC results for the first layer solid 
phase. Again, there is excellent agreement between the 
two methods. 

Following the typical procedure, we have determined 
the equations of state for the liquid and solid phases by 
fitting our energy values to the polynomials 

E/N = e Q + B{P^f + C{ P —^)\ (6) 
Pa Pa 

E/N = a + /3p + 7p 2 + Sp 3 . (7) 

for the liquid and solid phases, respectively. Values for 
the fitted parameters are given in Table |[ Errors for the 
liquid phase are indicated. The values shown for the solid 
phase were used in the plot of Fig. ^. The values given 
for /3,7, and 5 are accurate to one significant figure. 
For the liquid phase, the parameters eo and po are 
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FIG. 1. Energy values for the first layer liquid. The circles 
are ourrfesults obtained with PIMC. The squares are GFMC 
results.u The line is a least squares fit of the polynomial, Eq. 
([j]) to our data. 



-125.0 



-130.0 



-135.0 



-140.0 



Parameter Liquid 


Parameter Solid 


e (K) -141.689(14) 
po(A~ 2 ) 0.0450(06) 
B(K) 2.42(27) 
C(K) 3.29(67) 
X 2 /v 2.57 


a(K) -157.46 
f3{K/A 2 ) -679.83 
7(K/A 4 ) -10504.77 
r5(K/A 6 ) 61032.24 
X> 1.3 



-145.0 

0.065 0.075 0.085 0.095 0.105 0.115 
Coverage (A~ z ) 



TABLE I. Fitted parameters for the liquid and solid equa- 
tions of state. The errors in the last digits of the parameters 
for the liquid phase are given in parenthesis. 



the equilibrium energy per particle and the equilibrium 
coverage. The equilibrium density we obtain, 0.0450 
atom/ A 2 , is in reasonable agreement with the GFMC 
value of 0.0443. Below this density, the system enters 
gas- liquid coexistence in the thermodynamic limit. In 
simulation, phase separation will not occur immediately 
because of the finite cost of creating the phase boundary. 
The system will instead enter a "stretched" , or negative 
pressure, state. If one continues to decrease the cover- 
age, however, the phase separated state will eventually 
become the favored state, and the "stretched" state will 
"snap" , forming a droplet plus vacuum. For a simulation 
with a constant number of particles and variable area, 
this occurs when the derivative of the spreading pres- 
sure, P — p 2 (de(p)/dp), is zero. At the spinodal point, 
the sound velocity becomes imaginary and the compress- 
ibility diverges. From the equation of state, we determine 
that this occurs at 0.034 atom/A 2 . For comparison, the 
spinodal density has been calculated to be betwesa-Qi©31 
and 0.038 atom/A 2 for two-dimensional hclium.EiElt!3 

From the polynomial fits for the two layers, we can es- 
timate the regions of phase coexistence for the two solids 
by using the Maxwell double tangent construction. See 
Fig. ||. Since we have a constant number of particles 
and a variable density, the Maxwell construction is found 
from the common tangent of the energy per particle ver- 
sus atomic area (inverse coverage). We determine the 
liquid-solid coexistence region to be from 0.0675 to 0.0700 
atom/ A 2 . This is comparable .-to the ranges found for 
both two-dimensional heliumErLJ and the helium film,cJ 
although the onset of full solidification that we find is 
at a somewhat lower density. The solidification density 
determined in this manner is subject to error, since our 
fitted parameters are not determined with a high degree 
of precision. 



FIG. 2. Energy values for the first layer triangular solid. 
The symbols have the same meaning as in Fig. |^. The line is 
determined from the fit of Eq. wa) to our data. 



B. Other properties 

As we have seen, in the absence of corrugations, the 
first layer has liquid-gas, self-bound liquid, liquid-solid, 
and solid regions. The gas phase has zero density at zero 
temperature. We expect the liquid-gas phase to have two 
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FIG. 3. Liquid-solid coexistence regions determined with 
the Maxwell construction. The dashed-dotted and dashed 
lines are the solid and liquid equations of state, respectively. 
Circles with error bars are calculated energy values for the 
liquid. Squares with error bars show values calculated for 
the solid phase. The unbroken line is the coexistence line. 
The arrows indicate the beginning and end of the coexistence 
region. 
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FIG. 5. Static structure function for the liquid phase at the 
indicated densities, in atom/A 2 . The static structure function 
in the (01) direction at a typical solid density (top) is also 
shown. 
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FIG. 4. Radial distribution function for the liquid phase at 
the indicated densities, in atom/A 2 . 




FIG. 6. Contour plots of the probability density in the 
plane of the substrate. The densities shown are, top row, 
left to right, 0.0256 and 0.0421 atom/A 2 ; bottom row, left to 
right, 0.0689, and 0.0918 atom/A 2 . 
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FIG. 7. Density profiles near second layer promotion. The 
coverages, in atom/ A 2 , are indicated. Normalization is chosen 
so that integrating the profiles gives the coverage. 

regions, one in which phase boundaries between these 
phases form, and one in which the liquid phase is arti- 
ficially stretched. This is governed by the energy cost 
required to create the phase boundary. The metastable 
liquid phase is spatially indistinguishable from the equi- 
librium state, but the formation of droplets can be vi- 
sualized. In Fig. ^ we plot the radial distribution func- 
tion g(r) for the liquid at coverages in the droplet region, 
near equilibrium, and near the onset of liquid-solid co- 
existence. Here, r is the magnitude of the projection of 
the distance vector between two particles onto the plane 
of the substrate. In the droplet region, the long range 
tail drops below unity, indicating the system does not 
uniformly cover the substrate. Near equilibrium the first 
peak has about the same height as the droplet phase, 
but the tail goes to unity at long ranges. The peaks are 
located at the same value of r, indicating that the aver- 
age separation distance of nearest neighbor particles in 
the droplet is close to the equilibrium value. The high 
density liquid has noticeably more correlation. The peak 
height is 1.5, compared to 1.2 near equilibrium, and the 
peak's position has shifted from 4.2 to 3.6 A, as one would 
expect from compression caused by increasing the den- 
sity. These values are in agreement with- values reported 
for the helium liquid in two dimensions,!! 2 ] indicating that 
the first layer liquid is very two-dimensional in character. 

The results for the static structure factor S(k) for the 
liquid phase are shown in Fig. || and can be interpreted 
similarly. At the lowest coverage, S(k) swings upward for 
low values of k instead of going to zero. This indicates 
the presence of droplets and a nonzero compressibility. 
For the two uniform liquid coverages, S(k — > 0) — ► 0. 
The structure function is more peaked for the high den- 
sity liquid. Figure || also shows results for a typical solid 
coverage. The peak heights for the fluid near equilib- 



rium and for the dense fluid are 1.2 and 1.6, respectively. 
These values are in reasonable agreement with, but some- 
what below, the two-dimensional values. 

Finally, we show in Fig. || probability contours for the 
first layer for the various regions discussed above. The 
lowest coverage shown is below the spinodal point, and 
as expected incompletely covers the substrate. Near the 
equilibrium density, the substrate is uniformly covered. 
In the dense liquid, localization can be observed. This 
coverage is in the liquid-solid coexistence region. Finally 
the system enters a triangular solid phase at the highest 
density. 

Promotion to the second layer may be determined by 
examining the density profiles for the dense solid. These 
are shown in Fig. at the indicated coverages. The 
occupation of the second layer is clearly visible at the 
coverages 0.1200 and 0.1270 atom/A 2 , while the cover- 
age 0.1180 shows no evidence of promotion. Integrating 
the profiles for the two largest coverages to the mini- 
mum between the peaks (4.5 A) gives 0.116 and 0.119 
atom/A 2 -fnt-ihe. first layer coverage. Experimentally, 
estimatest§E3ElliJ of first layer completion range from 
0.112 to 0.120 atom/ A 2 so our value is. in good agree- 
ment. The recent GFMC calculations^ obtain a value 
between 0.115 and 0.118 for the beginning of second layer 
promotion. 



IV. RESULTS WITH CORRUGATED 
SUBSTRATE 

In this section we present results obtained by using 
a more realistic model of the graphite substrate. Most 
of the results shown in this section were obtained with a 
simulation cell with the dimensions 25.56Ax 22.14A. The 
commensurate density corresponds to 36 helium atoms. 
Periodic boundary conditions were used in the plane of 
the substrate. Finite-size effects were examined by re- 
peating some calculations using a 34.08A x 29.51A sim- 
ulation cell, for which the commensurate density corre- 
sponds to 64 particles. 

A. At And Below Commensurate Density 

The principal conclusion we obtain from these calcu- 
lations is that the low density, low temperature phase 
of the first layer consists of commensurate solid clusters, 
rather than a liquid phase. In an earlier publication^ we 
observed the commensurate phase and investigated its 
melting with both static structure and specific heat cal- 
culations. The calculated y3 x y3 commensurate solid 
phase and its melt are shown in Fig. [| At 3 K, the film 
has solidified into the commensurate structure. The solid 
forms a sub-lattice that contains one-third of the adsorp- 
tion sites. The remaining two-thirds of the sites form two 
equivalent sub-lattices that are unoccupied. Raising the 
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FIG. 8. Distribution plot at the commensurate density, 
0.0636 atom/A 2 , for T=2.99 K and T=4.0 K. The filled circles 
indicate graphite adsorption sites. 




FIG. 9. Probability distributions for 0.0566 atom/A 2 at 
T=1.0 (top) and 2.5 K (bottom). The filled circles give the 
locations of graphite potential minima. 



temperature to 4 K causes melting. At this temperature 
the helium atoms no longer inhabit a single sub-lattice 
of substrate adsorption sites. All adsorption sites will be 
occupied with equal probability if the simulation is run 
for a sufficiently long time. 

Next, we present direct evidence that the low den- 
sity (below the commensurate solid density of 0.0634 
atom/ A 2 ) first layer consists of a solid cluster surrounded 
by a low density vapor at low temperatures. No liquid 
phase forms and so there is no possibility for first layer 
superfluidity. These findings axe in contrast to the most 
recent proposal for this regionO. 

The presence of a solid with vacancies and phase sep- 
aration can be visualized with contour plots of the prob- 
ability distributions, shown in Fig. ^| for the indicated 
coverage and temperatures. At the lower temperature, 
1.0 K, the vacancies have condensed into a single bubble 
region, as can be seen in the top plot of Fig. g. Note that 
periodic boundary conditions are being used, so the va- 
cancy regions in the figure are connected. The vacancies 
move very slowly at this temperature; we found equilibra- 
tion times to be very long for the vacancies to condense 
if they were initially spread throughout the lattice. We 
thus calculated the energy for this system twice: first 



with the vacancies initially spread through the solid and 
then with the vacancies condensed. The condensed en- 
ergy was found to be lower. At higher temperatures, 
the vacancies acquire enough kinetic energy to leave the 
phase separated state and diffuse into the solid. As a 
result, vacancies can become isolated. This is illustrated 
at 2.5 K in the lower plot of Fig. ^. A series of probabil- 
ity distribution plots reveals that these vacancies move in 
the simulation, so the equilibration problem encountered 
at 1.0 K is not present at this temperature. We note that 
we still see evidence of phase separation in contour plots 
at 2.0 K, while heat capacity peaks from experiment seem 
to indicate a transition at 1.5 K. 

Figures [II] and [ll] further confirm that the low density 
region contains solid clusters at low temperatures and 
that these clusters exhibit the melting behavior discussed 
in the previous section. We have calculated distribution 
plots for densities as low as 0.0207 atom/A 2 and observe 
solid clusters at all densities. 

We have also attempted to place a vacancy in a solid 
cluster to see if the cluster could support an isolated va- 
cancy and at the same time be in equilibrium with the 
low density vapor. We found that at 0.0424 atom/A 2 and 
1.0 K, the vacancy was spontaneously expelled from the 
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FIG. 10. Distribution plot at 0.0530 atom/A 2 , T=2.5 K 
and T=3.0 K 

solid cluster during thermalization. We conclude that 
at low temperatures the solid clusters will not support 
isolated vacancies. 

Further evidence for solidification into the commensu- 
rate structure in the simulation comes from calculations 
of the static structure factor. Typical results for static 
structure factors for coverages at and immediately below 
the commensurate solid density have been reported in 
our earlier paperE 2 !. In addition, the melting of the com- 
mensurate solid phase can be determined from the tem- 
perature dependence of the static structure peak height. 
Melting is signaled by a significant drop in the peak 
height and a large statistical fluctuation in a peak value 
near the melting temperatures. The melting temperature 
determined in this manner are given in Table Q. The 
density dependence of the melting temperature is consis- 
tent with the experimental phase diagram, although our 
melting temperatures are slightly higher than the experi- 
mental values. Heat capacity measurements indicate that 
the commensurate solid melts at about 3 K, and the low 
density (< 0.045 atom/A 2 ) melting peaks are at about 
1.5 K. 

Another way of estimating melting temperatures is 
from the temperature dependence of the energy. This 
is shown in Fig. [l2] for various densities. These curves 




FIG. 11. Distribution plot at 0.0424 atom/A 2 , T=1.0 K 
and T=2.0 K 

possess inflection points that lead to specific heat max- 
ima when differentiated. These peaks indicate melting. 
Figure [l^ shows two sample calculations of the specific 
heat at the indicated densities. The peak height and 
location change dramatically with density. Melting oc- 
curs at about 1.5 K and 3.5 K for 0.0353 (filled cir- 
cles) and 0.0636 atom/A 2 (open squares) respectively. In 
Fig.[L3| the calculated specific heat for 0.0353 and 0.0636 
atoms/ A 2 is shown. For comparison we also present the 
experimental measure specific heat (solid lines) for den- 
sities of 0.0367 and 0.0634 atoms/AA 2 . The calculated 
specific heat has peaks somewhat above the experimen- 
tal values but are consistent with the melting behavior 
exhibited in the contour plots and in the static structure 
factor calculations. 



B. Maxwell Construction 

We determine phase ranges by applying the Maxwell 
construction to the total ground state energy. For a sys- 
tem at constant volume with a varying number of parti- 
cles, a region of phase separation will be signaled by an 
unphysical upward curvature in the total free energy's 
dependence on density. The upper and lower bounding 
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TABLE II. Estimates of melting temperatures from the 
temperature dependence of the static structure peaks. 



Coverage (atom/ A 2 ) 




0.0424 


2.0 


0.0530 


2.5 


0.0566 


3.0 


0.0636 


3.33 




-140.0 
-141.0 
-142.0 
-143.0 
z -144.0 

LLJ 

-145.0 
-146.0 
-147.0 
-148.0 

0.8 1.3 1.8 2.3 2.8 3.3 3.8 4.3 4.8 5.3 5.8 
Temperature (K) 

FIG. 12. Temperature dependence of the energy per par- 
ticle. The coverages are 0.0353 (circles), 0.0566 (triangles), 
and 0.0636 atom/A 2 (squares). 



densities of this region are connected by a common tan- 
gent line. The total free energy for all intermediate den- 
sities lies on or above this line, either because creating a 
surface between the two phases costs a finite amount of 
energy or because the system remains unphysically ho- 
mogeneous. In the thermodynamic limit the system will 
separate into the two phases at the bounding densities. 

The Maxwell construction at nonzero temperatures 
should be applied to the total free energy. This is not 
directly accessible from PIMC, however. Instead, we use 
a limiting process to determine effectively ground state 
energy values. All energy calculations are performed at 
low temperatures. The temperature is then raised and 
the energy is recalculated. If the two values are the same 
within error bars, we conclude that we have obtained ef- 
fectively zero temperature energy values. This allows us 
to apply the Maxwell construction to the total energy, 
since it is the same as the total free energy at zero tem- 
perature. 

We now can apply this procedure to the total energy 
calculations of the first layer solid. Total energy calcu- 
lations for a range of densities using a simulation cell 
designed to accommodate the commensurate solid are 



<D 4 



CD r> 
CO 








1 2 3 4 5 6 

Temperature (K) 

FIG. 13. Specific heat for 0.0353 (filled circles) and 0.0636 
atom/A 2 (squares). The dashed line is a guide to the eye. The 
solid lines are the experimentally determined specific heat at 
0.0367 and 0.0634 atom/ A 2 . 



TABLE III. Energy/particle versus coverage. The first col- 
umn gives the temperature of the calculation. The number in 
parenthesis gives the error in the last decimal place. 



T(K) 


a(A~ 2 ) 


E/N (K) 


1.00 


0.0353 


-143.73(8) 


1.00 


0.0424 


-144.03(9) 


1.00 


0.0530 


-144.31(7) 


1.00 


0.0566 


-144.56(7) 


1.00 


0.0636 


-145.12(8) 


1.33 


0.0707 


-142.71(8) 



shown in Fig. [Xj. The values for the energy per par- 
ticle are given in Table [II. The cell dimensions are 



25.560A x 22.136A, and the number of particles varies 
from 20 to 40. The dashed line in Fig. |l4| is a straight 
line connecting the lowest density accessible in our lat- 
tice and the density that corresponds to 1/3 coverage. 
Notice that all intermediate energy values are above this 
straight line. From the Maxwell construction, this in- 
dicates that the intermediate energy values are unstable 
and will phase separate into the two stable phases (the 
vapor and the commensurate solid) that bound the un- 
stable region. 

The binding energy for a single particle on the sub- 
strate can be easily calculated. We find this to be Eb = 
— 143 .09 ± 0.27 . This is comparable to the estimated 
values^ Eb = — 141. 75± 1.50 K from scatteringc^EJ and 
Eb = —142.33 ± 1.97 K from thermodynamic analysisEJ. 
Our value for the binding energy was calculated at 0.4 
K and confirmed to be the ground state value by the 
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-142 




0.02 0.04 0.06 0.08 

Density (atoms/A 2 ) 

FIG. 14. Total energy versus coverage. The dashed line is 
a straight line connecting the lowest density accessible in our 
lattice and the density that corresponds to 1/3 coverage. 

limiting procedure discussed above. By subtracting this 
energy from the energy per particle of the commensurate 
solid phase, we obtain the condensation energy per par- 
ticle for the two-dimensional solid, E^d = —2.03 ± 0.20 
K. This is comparable to, but slightly higher than, the 
two-dimensional energy (-1.06 K) for the commensurate 
phase found by the variational calculations ofc3 for the 
same interaction model. 



C. Above Commensurate Density 

Figures [II] and [l6] illustrate the domain wall solid and 
incommensurate triangular solid phases, respectively. 
The domain wall solid consists of patches of commensu- 
rate solid on different sub-lattices. Linear domain walls 
occur at the boundaries between these regions. In the 
incommensurate solid, the atoms form a triangular solid 
that does not have a periodic relationship with the un- 
derlying adsorption sites for lengths scales less than the 
minimum dimension of the simulation cell. 



V. TESTS OF OUR CONCLUSIONS 

The strength of our conclusions is limited by the ac- 
curacy of the interaction model that we use. First, it is 
possible that the substrate, way substantially alter the 
helium-helium interactiont^rca Inclusion of this effect, 
the so-called McLachlan interaction, has been shown to 
change the ground state phase from the commensurate 
solid to a low density liquid in two-dimensional vari- 
ational calculations. E3 We have repeated the low tem- 
perature density scans using the sam£_mediated interac- 
tions employed in these calculationsEJ. We found that 
while the energy per particle increases for all coverages, 




FIG. 15. Distribution plot of the domain wall solid at 
0.0742 atom/ A 2 (N=42 atoms) and T=1.0 K. Filled circles 
indicate adsorption sites. 

the commensurate solid remains the energetically favored 
phase. Second, another potential problem is that the 
helium-graphite interactions is too corrugated, thus fa- 
voring solidification. We have repeated the low temper- 
ature calculations at the commensurate density and at 
0.0424 atom/A 2 with the corrugation strength reduced 
by 25 %. The commensurate phase remains energeti- 
cally favored. Calculations were also performed with the 
corrugation reduced by 50 %, but it was found that the 
commensurate solid would not form for temperatures as 
low as 2 K, thus indicating the corrugations had been un- 
derestimated. We note finally that the melting behavior 
of the commensurate solid phase was not sensitive to the 
inclusion of the McLachlan term. 

Finally, we wish to discuss the arguments of Greywall 
and Busch (GB) against solid clusters and in favor of the 
superfluid phase. Their primary objection to solid-vapor 
coexistence is that this should be signaled by linear heat 
capacity isotherms for the entire region from zero cov- 
erage up to the commensurate density. Their published 
data shows that for temperatures from 0.2 K to 0.5 K, 
the isotherms are linear only between 0.025 and 0.060 
atom/A 2 . At 0.1 K, the upper endpoint is about 0.055 
atom/ A 2 . As a possible explanation, we suggest that the 
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FIG. 16. Distribution plot of the incommensurate solid 
phase at 0.0994 atom/A 2 (N=50 atoms) and T=2.0. Filled 
circles indicate adsorption sites. 



departure from linearity below 0.025 atom/A 2 is caused 
by the presence of multiple finite-sized clusters. At low 
densities, solid clusters nucleate around surface defects. 
Initially, there are many small metastable clusters with 
large perimeter-to-area ratios. Increasing the density in- 
creases the size of the clusters until the surface is covered 
by a few large solid clusters with negligible boundary ef- 
fects. Thus, the heat capacity exhibits linear behavior 
only after the solid clusters are sufficiently large so that 
the perimeter-to-area ratio is small. This presumably oc- 
curs for coverages above 0.025 atom/A 2 . GB have used 
a similar explanation in their arguments for solid-liquid 
and liquid-gas coexistence in regions that do not have 
linear isotherms. 

GB's identification of coverages near 0.04 atom/A 2 as 
liquid is based partly on simulation results for 2D helium 
on a flat substrate, the most relevant calculations then 
available. As GB note, the large peak associated with the 
melting of the uniform commensurate solid phase first 
emerges abcuse 0.04 atom/A 2 . 2D helium is a liquid near 
this densityE 2 ], suggesting that first layer coverages below 
0.04 may be liquid. Unlike the purely 2D simulations, our 
calculations take the role of substrate effects into account. 
As we have shown, surface corrugations push the density 



of the energy minimum up from about 0.04 on a flat 
substrate to 0.0636 atom/A 2 and produce solidification. 
GB also show that their low density heat capacity results 
are in general agreepicnt with a PIMC calculation for 
2D superfluid heliumc 2 ], suggesting that there might be a 
superfluid transition in the first layer. We have shown in 
Fig. [l3] that these rounded heat capacities are produced 
by the melting of a solid cluster and are not associated 
with a superfluid transition. 

In closing, we would like to remark-that a promising di- 
rection for monolayer superfluidityQeded lies with helium 
adsorbed on alkali substrates, particularly lithiumc3'E3. 
These substrates have much smaller corrugations and a 
much weaker attraction, allowing the first layer helium 
film to be a liquid. The phenomenon competing with 
superfluidity for these substrates is pre-wetting, rather 
than solidification. 



VI. SUMMARY 

Our first layer calculations were performed both with 
and without substrate corrugations. When corrugations 
are neglected, the first layer film resembles a purely two- 
dimensional film. We determined that the film consists of 
gas, liquid, and solid phases. These are separated by co- 
existence regions, and we determined the coverage ranges 
for all phases at low temperatures by using the Maxwell 
construction. The first layer liquid has an equilibrium 
density of 0.0450 atoms/A 2 . Below this density the sys- 
tem phase separates. This region is divided into an unsta- 
ble region, where liquid droplets form, and a metastable 
region in which the system over-expands instead of form- 
ing an interface. These two regions are separated by a 
spinodal point at 0.034 atoms/A 2 . At higher densities 
the system enters a narrow region of liquid-solid coexis- 
tence between 0.0675 and 0.0700 atoms/A 2 . Above these 
coverages the system is in a triangular solid phase. The 
beginning of layer promotion occurs between 0.116 and 
0.119 atoms/A 2 . All of our calculations are in agreement 
with recent Green's function Monte Carlo results.Ej 

When corrugations are included in the first layer, the 
phase diagram is substantially altered. We find that a 
v3 x \/3 commensurate solid occurs, in agreement with 
numerous experiments. By examining the temperature 
dependence of the static structure function and the spe- 
cific heat, we find that this solid melts at approximately 
3.5 K, compared with the experimental melting temper- 
ature of 3 K. We further find that the commensurate 
solid phase is energetically favored. At densities below 
commensuration, the system phase separates into com- 
mensurate solid clusters and a low density vapor. No 
liquid phase occurs at low temperatures. 
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